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In ideal MHD, shear Alfven eigenmodes may experience dissipationless damping due 
to resonant interaction with the shear Alfven continuum. This continuum damping 
can make a significant contribution to the overall growth/decay rate of shear Alfven 
eigenmodes, with consequent implications for fast ion transport. One method for 
calculating continuum damping is to solve the MHD eigenvalue problem over a suit¬ 
able contour in the complex plane, thereby satisfying the causality condition. Such 
an approach can be implemented in three-dimensional ideal MHD codes which use 
the Galerkin method. Analytic functions can be fitted to numerical data for equilib¬ 
rium quantities in order to determine the value of these quantities along the complex 
contour. This approach requires less resolution than the established technique of cal¬ 
culating damping as resistivity vanishes and is thus more computationally efficient. 
The complex contour method has been applied to the three-dimensional finite ele¬ 
ment ideal MHD code CKA. In this paper we discuss the application of the complex 
contour technique to calculate the continuum damping of global modes in tokamak 
as well as torsatron, W7-X and H-1NF stellarator cases. To the authors’ knowledge 
these stellarator calculations represent the first calculation of continuum damping for 
eigenmodes in fully three-dimensional equilibria. The continuum damping of global 
modes in W7-X and H-1NF stellarator configurations investigated is found to depend 
sensitively on coupling to numerous poloidal and toroidal harmonics. 

PACS numbers: 52.30.Cv, 52.35.Bj, 52.55.Fa, 52.55.He, 52.65.Kj 
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I. INTRODUCTION 


Shear Alfven modes are a well-established feature of toroidal magnetically confined plas¬ 
mas. These weakly damped modes are of substantial interest to fusion research as they may 
become unstable through resonant interaction with fast particles. This can degrade fast 
particle confinement, reducing plasma heating and damaging plasma facing components. 
Inhomogenous toroidal plasmas give rise to a continuum of modes of varying frequency, 
resonant on particular flux surfaces. Asymmetries in the magnetic geometry of the plasma 
lead to couplings between different poloidal and toroidal harmonics, resulting in gaps in the 
continuous spectrum in which global modes can occur. Examples of this type of gap mode 
include the toroidicity and elipticity induced shear Alfven eigenmodes (the TAE and EAE 
respectively) 1 .^. Global modes can also occur due to extrema of the continuous spectrum, 
which can arise from variation of the magnetic geometry and density of the plasma^. 

Resonant power absorption occurs where the eigenfrequency coincides with the frequency 
of a continuum branch. Physically, dissipative effects such as charge separation and mode 
conversion occur at these resonances, resulting in continuum damping. These effects and 
the resulting damping can be described using a kinetic theory of plasmas. The standard 
approach to calculate this damping numerically is to represent these processes using a resis¬ 
tive magnetohdrodynamic (MHD) model. In these models, continuum damping represents 
the limit of resistive damping as resistivity vanishes. However, this damping can also be 
found in non-dissipative ideal MHD theory. In this case the singularity resulting from the 
continuum resonance must be treated in accordance with the causality condition, which im¬ 
plies integration around the singularity to give a negative imaginary frequency component. 
This is formally analagous to the treatment of Landau damping in collisionless plasmas^. 
As they have complex frequencies, the solutions to the eigenvalue problem where continuum 
resonance interaction occurs are more accurately described as “quasi- modes”-. However, the 
more common description of the solutions as “modes” is used in this work. 

A number of techniques have been developed for calculating continuum damping using 
ideal MHD. Analytic treatments have been developed based on solution of the wave equa¬ 
tion using asymptotic matching^, representing the resonance as a perturbation following a 
ballooning transformation^ and representing the resonance as a perturbation to the La- 
grangian for a global model!. These techniques have been limited to the analysis of TAEs in 
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large aspect-ratio circular cross-section tokamaks. Continuum damping has previously been 
calculated numerically by including a set of analytically determined jump conditions due 
to singularities at continuum resonances^. Continuum damping has also been computed 
for shear Alfven eigenmodes in tokamak and stellarator cases using a numerical method 
in which the eigenvalue problem was solved along a complex contour—. However, in these 
cases, calculations were resticted to only two coupled harmonics in the large aspect ratio 
limit. Nevertheless, this complex contour technique proved more computationally efficient 
than the standard approach of finding damping in the limit of vanishing resistivity. 

In this paper we describe the results of a complex contour method for calculating contin¬ 
uum damping of shear Alfven eigenmodes in general three-dimensional toroidal geometry. 
Thus, modes involving couplings between many different poloidal and toroidal harmonics 
may be computed. This allows, for the first time, the capability to compute continuum 
damping in stellarators with realistic geometry. Section [XT] briefly describes the reduced 
shear Alfven wave equation and its solution using the code CKA—. The complex contour 
technique used to compute the continuum damping is discussed in Section IIIII This section 
outlines the choice of contour and calculation of equilibrium quantities for complex values 
of the radial coordinate. This technique is applied to a TAE in a tokamak with circular 
cross-section in Section IIV1 where it is verified against the previously developed resistive 
MHDM and large-aspect ratio complex contour— continuum damping calculations. The 
complex contour technique is used to calculate the continuum damping of TAEs and EAEs 
in stellarators with torsatron, liclias and heliac configurations in Section [V] 


II. SHEAR ALFVEN EIGENMODES 

CKA (Code for Kinetic Alfven waves) is a finite element MHD code which solves reduced 
MHD equations for <f>. the peturbation to the scalar potential—. This code can be used to 
find Alfven eigenmodes in an incompressible plasma in either two or three dimensions. The 
code uses Boozer coordinates (s,9,(j)), where s e [0,1] is the square root of the normalised 
toroidal flux s = >/ . CKA represents a discretisation of the following Alfven wave 
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equation using the Galerkin method 
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Here Va = ^^ Uiin represents the Alfven speed, where n* is the number density of ions and 
ro ( : is their mass. Pressure is represented by P(°), magnetic held strength by B and held 
line direction by the unit vector b. The component of the current parallel to the held lines 
is j ,|°\ In this paper negligible /3 is assumed, so that terms involving P(°) in equation (fTi) 


can be neglected. Low ion and electron temperatures are considered so that pi = ^ k '^ Ti 
and p s = v/A 'gg Te are much smaller than the scale length of the waves Rq. Therefore the 
second term on the left hand side of equation (JT]) which represents the parallel electric held 
component, can be ignored. Frequency can also be expressed as a dimensionless quantity 

normalised to the Alfven frequency, fl = — where Rq is the mean major radius of the 

v A Ro 

plasma and is the Alfven speed at the magnetic axis. This equation represents an 
eigenvalue equation of the form 


L 1 [<h(s)]=u; 2 L 2 [$( S )]. 


where L\ and L 2 are self-adjoint differential operators for real s. 


( 2 ) 


Solutions to the eigenvalue equation are represented using B-splines, B i} to represent 
variation in each dimension^. B-splines are splines with support over a minimal number of 
knot intervals for a given order. Thus for a two-dimensional (axisymmetric) case, the scalar 
potential perturbations are represented as 

N s N e 

4'M) = EZ<VB,M BjW, (3) 

i=l j=1 
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while in the three-dimensional case, these perturbations are represented as 

N s Nq N<f) 

(M,0) = E E E c ljk B t (s) B j ( e ) B k (0). (4) 

i= 1 j =1 fc=l 

These splines are defined such that knots lie on a rectangular gridiS. Using the Galerkin 
method, the eigenvalue problem can be discretised by first substituting the above into 
equation ([I]). The resulting expression is then multiplied by Bi (s) Bj ( 9 ) (two-dimensional 
case) or Bi (s) Bj ( 6 ) B k (0) (three-dimensional case) and integrated by volume. Assigning 
a unique index l to each set of i, j and k the resulting equation can be expressed as a matrix 
eigenvalue problem. The numerical solution represents a projection of the exact solution 
<f> onto the space spanned by the basis functions. 

Eigenmodes are clearly periodic in the poloidal and toroidal directions. Moreover, they 
are often dominated by a small number of harmonics with close poloidal and toroidal mode 
numbers. Therefore it can be advatageous to use a phase factor extraction method. In this 
case the numerical solution is expressed as e }( me + n <t>) where rh and n are poloidal 

and toroidal mode numbers and Nf p is the number of toroidal field periods. The eigenvalue 
problem is reformulated in terms of which is a linear combination of B-splines. This 
reduces the poloidal and toroidal resolution required for the spline basis. 


III. COMPLEX CONTOUR METHOD 

Continuum resonance damping can be calculated by solving the TAE wave equation 
along a complex contour in the radial coordinate^. This is analogous to the analysis of 
Landau damping of plasma oscillations^. In that case dissipationless damping arising from 
inhomogeneity in velocity space is analysed by choosing an appropriate contour in velocity 
space. Poles due to continuum resonances are circumvented, and therefore singularities in the 
solution are avoided. Damping is represented by an imaginary component of eigenfrequency. 
In order to implement this technique the equilibrium quantities in equation (JTD must be 
analytically continued to determine their values at complex values of the radial coordinate 
along the integration contour. The damping of TAEs computed using the complex contour 
technique in ideal MHD has previously been shown to converge to the same value as in the 
limit as resistivity is reduced to zero in resistive MIID-ii. 
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The causality condition specifies that poles due to continuum resonances should lie be¬ 
tween the integration contour and the real axis. Physically, this condition ultimately arises 
from the requirement that any perturbation to the state of the plasma preceedes its effect 
on other plasma parameters. This in turn implies that modes decay rather than grow with 
time and hence that 3 (u) < 0. The imaginary frequency component causes the poles to 
move from real to complex values of the radial coordinate. 

If the integration contour is shifted away from poles due to continuum resonances, the 
calculated modes will vary less rapidly in their vicinity. Thus, the solution to the eigenvalue 
problem can be accurately represented by basis functions defined on a coarser radial grid. 
Consequently, the eigenvalue will have faster convergence with respect to radial grid reso¬ 
lution for greater deformations of the integration contour. However, excessive deformations 
may result in the disappearence of Alfven eigenmodes as the continuum gaps they exist 
in close. In such cases the chosen complex contour intersects with analytically continued 
branches of the continuous spectrum. Additionally, non-physically meaningful poles may 
result due to the analytic continuation of equilibrium quantities. These poles should not lie 
between the integration contour and real axis, in order to avoid spurious contributions to 
the eigenfrequency. 

The continuum damping calculations described below are performed using a complex 
contour described by: 

s - x + if (x). (5) 

Here, x e [0,1] is a real parameter. The real valued function / (x) is either chosen to be a 
quadratic: 

/ (x) = ax (x - 1) (6) 

or the product of a quadratic and a Gaussian function 

/ X—X'y G 

/ (x) = ax (x - 1) e V x p ) . (7) 

The parameters a, xp and x 7 respectively determine the scale, width and location of the 
deformation into the complex plane. Defining a profile using the latter function allows a more 
localised contour deformation than the former. Thus, larger deformations in the vicinity of 
the pole can occur while avoiding non-physically meaningful poles. Therefore, it is possible 
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to calculate continuum damping with a smaller grid resolution. However, as the deformation 
becomes narrower and its scale increases, the rate of variation of the solution in that region 
is expected to increase. Thus, increased radial resolution may be required where a is too 
large and xp is too small. 


A. Equilibrium quantity representation 

In order to solve the eigenvalue problem in equation (JT]) it is necessary to calculate values 
of equilibrium quantities for complex values of s. These include density p, toroidal flux ijj, 
poloidal flux, fj p , current j, magnetic field, B, pressure P(°), Jacobian yjg and metric tensor 
elements g l] . One possible approach to estimating these quantities for complex s is by using 
a truncated Taylor series expansion. Where the radial coordinate s has real component s r 
and imaginary component s*, expanding the function y ( s ) about s r yields 

y(s) = y(s r ) + if (s r ) y' (s r ) - ^ (/ (s r )) 2 y" (s r ) - ^ (/ (s r )) 3 y'" (s r ) (8) 

This approximation may be used with spline representations of equilibrium quantities. How¬ 
ever, these piecewise polynomial functions are not analytic, which can lead to branch cuts 
when extended to the complex plane. If the complex contour deviates sufficiently from the 
real axis, the calculated values of the functions might contain significant contour dependent 
errors. Integration across branch cuts may thus affect the result of the finite element calcu¬ 
lation, leading to poor eigenvalue convergence with respect to the chosen complex contour. 
This is found to be the case for both the axisymmetric tokamak and stellarator examples, 
when the first three terms in the expansion for each quantity were retained. In these cases 
poor convergence is observed with respect to the choice of contour (as parameterised a, xp 
and x 7 for contours described by equation (171) ). 

Consequently, in these calculations real and imaginary components of equilibrium quan¬ 
tities along the chosen complex contour are each represented in terms of cubic B-splines. 
These piecewise representations are determined from analytic representations of the equi- 
libium quantities. Equilibrium data from the VMEC code^ are used to obtain a Fourier 
series in 6 and f of polynomial functions in terms of s. These polynomial functions are calcu¬ 
lated based on the equilibrium data using a least squares fitting routine, ignoring unreliable 
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data close to the magnetic axis. The use of polynomial functions to represent equilibrium 
quantities allows trivial analytic continuation from the real axis to determine their values at 
complex s. The aforementioned spline representations of the non-flux function equilibrium 
quantities are calculated by evaluating these series at positions along the chosen complex 
contour in s and at various values of 9 and <f>. Use of spline representations for these quantities 
is advantageous as the evaluation of polynomials is relatively computationally inexpensive. 


IV. TOKAMAK RESULTS 

The complex contour technique for calculating continuum damping was applied to a TAE 
in an axisymmetric tokamak with circular cross-section. An inverse aspect ratio of e = 0.225 
is chosen for this case. The rotational transform is described by the polynomial profile 
l 0) = 0.95016 - 0.67944s 2 + 0.62286s 4 - 0.41244s 6 + 0.1219s 8 + 0.0042185s 10 - 0.0013979s 12 . 
Here the rotational transform is defined as i = liniA^oo where A 9 and A (ft are respectively 
the displacement in the poloidal and toroidal direction along a magnetic held line. The 
density profile is chosen to be 

^i(s) = -y(l-tanh^^-^JJ, (9) 

where Ai = 0.8 and A 2 = 0.1. This approximates the choice of t and rti used by Konies and 
Kleiber—. The functions n* (s) and t (s) are overlaid on the resulting continuous spectrum 
shown in figure CD This equilibrium results in a TAE due primarily to the coupling of the 
(m,n) = (2,-2) and (3,-2) harmonics. Complex contours of the form indicated in equation 
(J6]) were used. 

The frequency of this TAE was calculated for various values of the contour deformation 
parameter a and radial grid resolution N s . For this axisymmetric plasma the solution 
is represented by splines defined on a two-dimensional rectangular grid in s and 9. The 
results of this calculation are plotted in figure [21 demonstrating convergence of continuum 
damping with respect to these parameters. Thus, a complex normalised TAE frequency of 

= 0.39515 - 0.00174465z is obtained, resulting in a damping ratio of = -0.00441515. 
This convergence indicates that sufficiently large values of a and N are used that the mode 
structure is resolved. The most demanding element of this structure in terms of radial grid 




FIG. 1. Shear Alfven continuous spectrum as a function of radial position for a circular-cross section 
tokamak computed using CKA. The avoided crossing between the (m = 2 ,n = -2) branch (o) and 
the (3,-2) branch (o) result in a TAE (not shown). The rotational transform profile, i, (dashed 
line) and normalised density profile, (solid line) are also plotted with the scale shown on the 
right axis. 


resolution is the continuum resonance interaction, particularly where the chosen contour 
passes close to the pole due to the continuum resonance. In that region there is rapid 
variation in the perturbation to the scalar potential due to the mode, requiring a large 
density of grid points to resolve. A reduction in the poloidal grid resolution from Ng = 20 to 
10 results in a change in damping ratio of less than 0.07%, indicating satisfactory convergence 
with respect to this parameter also. 
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FIG. 2. Convergence of continuum damping with respect to the contour deformation parameter a. 
Damping ratio is plotted as a function of a for N s = 400 and Ng = 20 (■), N s = 800 and Ng = 20 
(•), N s = 1200 and Ng = 20 (▼) and N s = 1200 and Ng = 10 (A). 


Continuum damping of Alfven eigenmodes in axisymmetric tokamaks has been com¬ 
puted previously by determining the limit of resistive damping as resistivity tends to zercM. 
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Therefore, the resistive calculation can be used as a benchmark for the complex contour 
calculation. This calculation is implemented in CKA by inserting an additional term on the 
left hand side of the equation ([[]), 

($(*))) ( 10 ) 

Damping is plotted as a function of radial grid resolution for a series of different values of the 
artificial damping parameter 5 in figure [3j This plot indicates convergence in the damping 
ratio with respect to both 5 and N. Complex normalised TAE frequency is estimated to be 
D = 0.395151-0.0017452D corresponding to a damping ratio of = -0.00441658. In this case 
the density of the grid in the radial coordinate is increased in the vicinity of the continuum 
resonance at s = 0.964, to increase the rate of convergence. Satisfactory convergence can 
also be demonstrated with respect to the poloidal coordinate, as the change in the damping 
ratio computed when the poloidal grid resolution is reduced from Ng = 20 to 10 is less than 
0.07%. 

The small difference between the damping calculated using the complex contour and 
artificial resistivity methods can be attributed to the errors present in each. Approximation 
of equilibrium quantities in terms of cubic splines may account for the variation in the 
continuum damping ratio with contour deformation a. In the artificial resistivity method a 
finite resistivity is required to resolve the resonant interaction with the continuum due to the 
finite radial grid resolution^. This results in a small resistive damping occuring throughout 
the plasma in addition to the continuum damping very close to the magnetic surface on which 
the continuum resonance occurs. Reducing the resistivity to improve accuracy requires a 
greater radial grid resolution to obtain convergence, as the finite width of the peak in the 
mode due to continuum resonance must be adaquately resolved. For sufficiently large a the 
damping found using the complex contour technique is found to converge more rapidly with 
respect to N s than for the resistive technique. 

It is also interesting to compare the continuum damping calculated in the above fully 
two-dimensional case to that found using the large aspect-ratio approximation. This ap¬ 
proximation results in a differential equation derived by Berk et al. to describe TAEs^. 
The application of the complex contour technique based on this formula has been described 
previously by Konies and Kleiber— and Bowden et alM. Using the shooting method de- 
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FIG. 3. Convergence of damping as the artificial resistivity parameter 6 is reduced for a circular 
cross-section tokamak. Damping ratio is plotted as a function of N s for 6 = 5 x 10 -9 (■), 1 x 10 -8 
(•), 5 x 10 -8 (▼) and 1 x 10 -7 (A). 


scribed in the latter and considering only the (2, -2) and (3, -2) harmonics, a TAE is found 
with normalised frequency hi = 0.370887 - 0.00174204* corresponding to a damping ratio 
of = -0.00469695. Thus the calculation based on the large aspect ratio approximation 
agrees reasonably closely with the fully two dimensional calculation in this case (to within 
6.6%). This suggests that modifications to TAEs due to terms of higher order with respect 
to aspect ratio and couplings to additional harmonics do not have a large effect on their 
continuum damping in this case. 

The three-dimensional continuum damping calculation was verified through application 
of the fully three-dimensional calculation to this tokamak case. That is, the quadratic form 
was computed using integrals for basis splines defined on a three-dimensional grid. In the 
resulting solution, toroidal variation of eigenmodes is represented by splines with a radial 
grid resolution of = 10 in this dimension and N s = 400 and Ng = 10. This results in a TAE 
frequency which agrees with the corresponding two-dimensional calculation to 7 significant 
figures, the precision provided in the output of the code. 


V. STELLARATOR RESULTS 
A. Torsatron eigenmodes 

The continuum damping of a TAE and an EAE in a torsatron configuration were calcu¬ 
lated. The rotational transform for this case is defined by the polynomial t (s) = 0.4319 + 
0.23407s 2 + 0.042125s 4 + 0.008341s 6 . This profile is typical of torsatron configurations, with a 
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minimum near the magnetic axis and strong magnetic shear near the edge!!. For simplicity, 
the same density profile is used here as in the tokamak case in Section IIVI The stellarator 
has Nf p - 20 toroidal held periods with a major radius of R = 20m. The plasma cross section 
is elliptical, with semi-major axis of a = 1.13m and semi-minor axis of b - 0.867m. TAEs 
have previously been observed in torsatrons such as the Large Helical Device (LHD )!& and 
the Compact Helical System (CHS)^^. In some instances continuum damping has been 
identified as the dominant cause of damping for these modest. In the present case a large 
EAE gap exists due to the elongation of the plasma, while the TAE gap is small due to its 
large aspect ratio. The continuous spectrum is plotted as a function of radial position in 
figure m demonstrating these features. 



FIG. 4. Shear Alfven continuous spectrum as a function of radial position for the torsatron case 
computed using CKA. The avoided crossing between the (m = 3 ,n - -2) branch (o) and the (4, -2) 
branch (o) produce a TAE in the resulting spectral gap. Likewise the avoided crossing between 
the former branch and the (5, -2) branch (o) gives rise to an EAE in the corresponding spectral 
gap. The locations of the maxima of these global modes are shown, approximately corresponding 
to the radial locations of the gaps. The rotational transform profile, l, (dashed line) and normalised 
density profile, 4H-, (solid line) are plotted with the scale shown on the right axis. 


A TAE was found due to coupling between the (3,-2) and (4,-2) harmonics. This TAE 
has resonances with n - 2 branches of the continuum at r & 0.28 and r & 0.94. Using 
complex contours described by equation (J7J) it is possible to calculate the contribution of 
each resonance to the overall damping separately by localising the contour deformation near 
that resonance. This approach assumes that the other continuum resonance interaction does 
not significantly alter the mode structure. This is not always the case, even for low values 
of damping!!!. However, in this example we will assume that damping due to the resonance 
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near the edge is much larger than that near the core. This is expected, as the avoided 
crossing of the (3,-2) and (4,-2) continuum branches is closer to the resonance with the 
(4,-2) branch, so the mode would be expected to have greater amplitude at this crossing. 
Moreover, the amplitude of the (4, -2) harmonic is expected to be far larger than the (5, -2) 
harmonic throughout the plasma. Thus, = 0.96 is initially chosen, localising the contour 
deformation near the resonance close to the edge of the plasma. 
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FIG. 5. Convergence of continuum damping ratio for a TAE due to coupling of the (3,-2) and 
(4, -2) shear Alfven wave harmonics in a torsatron plasma, with respect to (a) radial grid resolution 
N s , (b) poloidal grid resolution Ng and (c) toroidal grid resolution AG. In each case, two of N s = 300, 
Ng = 20 and AG = 8 are fixed and the other parameter is varied. Complex contours are chosen with 
parameters xp ~ 0.1 and x 7 = 0.96. Damping ratio is plotted for a = 0.1 (■), a = 0.2 (•), a = 0.5 
(▼) and a = 1.0 (A). 

The damping ratio of the TAE is plotted as a function of N s , Ng and for various 
values of a, xp = 0.1 and x 7 = 0.96 in figure [5] The ratio is found to be approximately 
constant with respect to each of these parameters, indicating convergence with respect to 
grid resolution in each dimension and the complex contour chosen. The most converged 
estimate of the complex normalised TAE frequency obtained is Q - 0.29031 - 0.00333875* 
corresponding to a damping ratio of = -0.0115006. When the complex deformation is 
localised near the inner resonance and grid resolution is kept constant, comparatively weak 
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damping is found. For a - 0.1 to 0.5, xp = 0.05 and x 7 = 0.28 a complex normalised 
frequency of 12 « 0.291398-2.7 x 10 -8 * and a damping ratio of » -9.2 x 10 -8 are computed. 
This validates the assumption that the continuum interaction with the (5,-2) branch is not 
significant for this mode. 
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FIG. 6. Convergence of continuum damping ratio for an EAE due to coupling of the (3,-2) and 
(5, -2) shear Alfven wave harmonics in a torsatron plasma, with respect to (a) radial grid resolution 
N s , (b) poloidal grid resolution Ng and (c) toroidal grid resolution Np. In each case, two of N s = 300, 
Ng = 20 and Np = 8 are fixed and the other parameter is varied. Complex contours are chosen with 
parameters xp =0.1 and x 7 = 0.96. Damping ratio is plotted for a = 0.1 (■), a = 0.2 (•), a = 0.5 
(▼) and a = 1.0 (A). 


The EAE for which continuum damping is calculated is due to coupling of the (3,-2) 
and (5,-2) harmonics. This mode has resonances with the (3,-2) continuum branch at 
s & 0.78 and with the (4, -2) continuum branch at s & 0.96. Once again, the contribution 
of each continuum resonance to the damping can be computed separately. As expected, 
it is interaction with the latter branch that contributes most to continuum damping. The 
convergence of the damping ratio due to the outer resonance, obtained using x 7 = 0.96 is 
plotted in figure [6l The converged case corresponds to a complex normalised frequency of 
12 = 0.48197 - 0.00248407* and hence a damping ratio of = -0.00515400 for the outer 
resonance. In contrast, examining the effect of the inner continuum resonance by setting 
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x 7 = 0.78 and varying between a - 0.1 and a = 0.5 results in a complex normalised frequency 
of Q = 0.484f7f - 0.0000149544i and a damping ratio of just ~ -0.000031. Thus, the 
contribution of the inner continuum resonance to the damping of the mode is very small 
compared with the outer one and can be safely neglected. 

The rapid convergence with respect to Ng and N,p in both the TAE and EAE cases 
indicates that there is relatively little coupling to poloidal and toroidal harmonics beyond 
the main two indicated indicated for those modes. This is confirmed based on the mode 
structures obtained in each case. The largest harmonics in the Fourier decomposition of the 
TAE and EAE are shown in figure [7] and figure 0 respectively, ft is found that the mode 
structure obtained on the complex contour is similar to that obtained solving the problem 
for real s. The largest differences in the amplitude and phase of the modes calculated in 
each case are observed for the (3, -2) harmonic. These differences can be seen to extend 
well beyond the region where the contour deformation is localised, and are thought to be an 
effect of the continuum resonances for that harmonic. 



FIG. 7. Largest Fourier harmonics for the TAE identified in the torsatron case as a function of 
contour parameter x. The dominant (3, -2) (blue) and (4, -2) (green) harmonics are labelled. 
The mode calculated along a complex contour of the form indicated in equation [7] with x a = 0.1, 
xp = 0.1 and x 7 = 0.96 (thick line) is compared with that calculated for real s (thin line). The real 
component is indicated by the solid line and the imaginary component is indicated by the dashed 
line. 
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FIG. 8. Largest Fourier harmonics for the EAE identified in the torsatron case as a function of 
contour parameter x. The dominant (3, -2) (blue) and (5, -2) (green) harmonics are labelled. 
The mode calculated along a complex contour of the form indicated in equation [7] with x a = 0.1, 
xp = 0.1 and = 0.96 (thick line) is compared with that calculated for real s (thin line). The real 
component is indicated by the solid line and the imaginary component is indicated by the dashed 
line. 

B. Helias eigenmode 

Similarly, a TAE was identified and its continuum damping calculated for a W7-X stel- 
larator configuration. This device is an Nf p - 5 toroidal field period helias which has a 
relatively complicated magnetic geometry produced by modular coils. Consequently, many 
different Fourier components contribute significantly to equilibrium quantities. These cause 
significant couplings between numerous different harmonics. However, as the rotational 
transform i remains close to unity throughout the plasma, only avoided crossings due to 
continuum branches with relatively high (m, n) occur in W7-X. While these modes involve 
larger poloidal and toroidal harmonics, the relatively low magnetic shear allows such modes 
to have a significant radial extent. 

A TAE is found which is primarily due to coupled (18, -16) and (19, -16) harmonics. The 
global mode calculation is done for a density profile defined using the equation in the Sec¬ 
tion 033 with parameters changed to = 0.4 and A 2 = 0.2. This density profile is expected 
to reduce the effect of poloidal mode couplings compared with that used in the tokamak and 
torsatron cases, as the continuum resonance is moved inward. Thus the required poloidal 
resolution is reduced, making the calculation more tractable. By reducing Ai the resonance 
is shifted inward to where the TAE amplitude is greater and by increasing A 2 the radial vari- 
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FIG. 9. Shear Alfven continuous spectrum as a function of radial position for the W7-X case 
computed using CKA. The avoided crossing between the (m=18,n = -16) branch (o) and the 
(19, -16) branch (o) leads to a TAE in the resulting spectral gap. The locations of the maxima of 
these global modes are shown, approximately corresponding to the radial locations of the gaps. 

ation in continuum frequency is decreased. A TAE mode and the corresponding branches 
of the continuous spectrum are shown in figure [UJ The TAE has a continuum resonance 
with the (18,-16) branch at s » 0.67. Although other continuum resonances are expected 
near the edge of the plasma due to the upward shift in continuum frequencies resulting 
from decreasing plasma density, it is believed that this resonance represents the dominant 
contribution to the continuum damping of this mode. This assumption is supported by the 
continuum damping calculations for the torsatron, where it is shown that continuum damp¬ 
ing is dominated by interaction with branches corresponding to the dominant harmonics of 
the global mode. Other TAEs with the same dominant harmonics and a different number 
of nodes are found in the gap. 

The computed damping ratio for the selected TAE is plotted as a function of N s , Ng and 
in figure [TO] These plots show convergence with respect to the number of splines used 
in each dimension, as well as the contour deformation parameter a. Thus the normalised 
complex frequency is estimated to be Q = 0.394342-0.00047676i based on the value obtained 
for the highest grid resolution in each dimension and smallest contour deformation. This 
corresponds to a continuum damping ratio of J- = -0.0012090001. Fine poloidal and toroidal 
resolution is required in this case due to strong couplings between the ( m,n ) and ( m',n') 
harmonics for a range of different m' - m and n' - n. This is in contrast to the torsatron 
case where couplings for m' -m - ±1 and n' -n = 0 harmonics (toroidicity induced coupling) 
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FIG. 10. Convergence of continuum damping ratio for a TAE due to coupling of the (18,-16) and 
(19, -16) shear Alfven wave harmonics in a W7-X plasma, with respect to (a) radial grid resolution 
N s , (b) poloidal grid resolution Ng and (c) toroidal grid resolution AC. In each case, two of N s = 150, 
Ng = 40 and Nq = 12 are fixed and the other parameter is varied. Complex contours are chosen 
with parameters xp = 0.1 and = 0.65. Damping ratio is plotted for a = 0.2 (■), a = 0.5 (•) and 
a = 1.0 (▼). 


and m' - m - ±2 and n' - n - 0 harmonics (elipticity induced coupling) are expected to be 
most significant and damping can be computed relatively accurately for lower poloidal and 
toroidal resolution. Thus we conclude that the continuum damping in the W7-X case is 
significantly affected by coupling to a range of non-dominant harmonics. 


C. Heliac eigenmode 

Finally, continuum damping was calculated for a non-conventional global Alfven eigen¬ 
mode (NGAE)2! in an H-1NF stellarator configuration^. This device is an Nf p = 3 toroidal 
held period flexible heliac. The magnetic geometry of H-1NF can be adjusted based on the 
ratio of currents through poloidal and toroidal held coils Kh- Here we consider a case for 
which Kh = 0.33. Again, the density prohle used for the tokamak case described in Section [TV] 
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is used. The continuous spectrum obtained for the H-1NF case is shown in figure fill 



FIG. 11. Shear Alfven continuous spectrum as a function of radial position for the H-1NF case 
computed using CKA. The maximum of the (m=10,n = -13) branch (o) near the edge of the 
plasma results in a NGAE, which has its maximum at the location indicated. The (8,-10) branch 
(o) is also shown. The rotational transform profile, i, (dashed line) and normalised density profile, 
(solid line) are plotted with the scale shown on the right axis. 


An NGAE is found which predominantly comprises the (10,-13) harmonic. This type 
of mode is associated with a local maximum of the corresponding branch of the continuous 
spectrum. In this case the mode is localised near the edge of the plasma, where density 
declines sharply and i approaches the rational surface. The NGAE is found to be resonant 
with the (8,-10) branch of the continuum at s rs 0.81. The convergence of the continuum 
damping due to this resonance is shown in figure [12] It is found that the normalised complex 
frequency of the mode is 12 = 0.277665 - 0.00162422z, corresponding to a damping ratio of 
■2- = 0.006647862. As for the W7-X TAE case, a relatively large poloidal and toroidal 
resolution is required due to the range of couplings induced by the complicated magnetic 
geometry. 


VI. CONCLUSION 

Continuum damping of Alfven eigenmodes in complicated three-dimensional geometries 
can be calculated by solving the ideal MHD eigenvalue problem along a complex contour 
in s. This is achieved by accurately fitting analytic functions to numerical equilibrium 
data and analytically continuing these to find values for complex s. An appropriate choice 
of complex contour allows relatively rapid convergence of this damping with respect to 
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FIG. 12. Convergence of continuum damping ratio for an NGAE due to the (10,-13) shear Alfven 
wave harmonic in an H-1NF plasma, with respect to (a) radial grid resolution N s , (b) poloidal grid 
resolution Ng and (c) toroidal grid resolution Np. In each case, two of N s = 96, Ng = 48 and Np = 16 
are fixed and the other parameter is varied. Complex contours are chosen with parameters xp = 0.1 
and = 0.8. Damping ratio is plotted for a = 0.1 (■), a = 0.2 (•) and a = 0.5 (▼). 


radial grid resolution. Coupling to numerous poloidal and toroidal harmonics can affect 
the damping of gap modes such as TAEs and EAEs in stellarators, as demonstrated by the 
convergence of damping with respect to grid resolution in the poloidal and toroidal direction. 

Analytic continuation of the eigenfunctions computed along complex contours could be 
used to determine the structure of modes with continuum resonances for real values of 
s. These resonances are associated with sharp peaks in the amplitude of <h and sudden 
changes in its phase. Conventional ideal MHD codes find a sharp peak at the resonance 
but do not reproduce the jump condition at this singularity which can be found by complex 
contour integration. Properly incorporating the effect of continuum resonances on the mode 
structure would enable more accurate computation of fast particle drive and transport by 
incorporating the potentially significant effect of continuum resonances. 
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